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We perform numerical simulations of the critical gravitational collapse of a spherically symmetric 
scalar field in 6 dimensions. The critical solution has discrete self-similarity. We find the critical 
exponent 7 and the self-similarity period A. 
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I. INTRODUCTION 

Critical behavior in gravitational collapse, as first found by Choptuik ||^, occurs at and near the threshold of black 
hole formation For a one-parameter family of initial data slightly above the threshold, the mass of the black hole, 
Mbh, scales like (p — p*y . Here p is the parameter, is its critical value and 7 is a constant that depends on the 
type of matter, but not on the family of data. For initial data slightly below the threshold, the maximum curvature 



• scales like {p * — p)~ where 7 is the same constant as in the black hole mass scaling law |^]. (The definition of the 
' scaling exponent 7 must be slightly generalized in n spacetime dimensions, where Mbh has dimension (length)""'^; 

• see section II.) The critical solution (p — p*) has either continuous self similarity or discrete self similarity, depending 
I on the type of matter. 

0^ , While critical gravitational collapse has been studied in many types of matter, the work has, in general, been done 

' in 4 spacetime dimensions. (The exception is work on analogs of gravitational collapse in 3 spacetime dimensions [Q). 

O One might therefore wonder whether critical behavior occurs in gravitational collapse in n spacetime dimensions for 

n > 4, and if so, how the properties of the critical behavior depend on n. 

^ , In this work, we perform numerical simulations of the collapse of a spherically symmetric scalar field in 6 spacetime 
dimensions. We find that the critical solution has discrete self-similarity. We find the scaling exponent 7 and the 

^ ' self-similarity period A. Section 2 reviews the Schwarzschild solution in n dimensions and shows how the definition 

• r of the 7 exponent should be generalized to the n-dimensional case. Section 3 gives the equations for the evolution of 

' the scalar field and the metric in a form suitable for our numerical simulations. The numerical method is presented 

^ in section 4 and our results in section 5. 



II. SCHWARZSCHILD IN N DIMENSIONS 

The Schwarzschild metric in n spacetime dimensions is |^ 

ds' = —)dt' + {l-- — —) dr^+r^dni_^ 1 

(n - 2)A„_2 f ("'-2)^1-2?' 

where ^„-2 = 27r^"~-^)/^/r(^i^) is the area of the unit {n — 2)-sphere. The quantity M in the metric (||) is the 
spacetime mass. M is normalized so that when dll) is the exterior metric of a nearly Newtonian static fiuid ball. 
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then M is J ToodV. With this normahzation, the force on a static unit-mass particle at distance r approaches 

(|^i^M/r"-4 as r ^ oc. 

It is important to note that M has dimension (length)" ^. This implies the following scaling behaviors for critical 
phenomena in n dimensions. If below threshold the maximum curvature [which has dimension (length)"^] scales like 
(p * — then above threshold Mbh scales as (p — p*)'^""^^'''. 



III. EVOLUTION EQUATIONS 



Our matter model is the one used in reference Q: a massless, minimally coupled, self gravitating scalar field 
The stress-energy of the field is 



Tab = ^a<t>^b 
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(2) 



In n spacetime dimensions, a spherically symmetric metric has the form 



ds^ = -e^''du^ - 2e''+^dudr + r^d^^ 



(3) 



where dT,^ is the metric of the unit (n — 2)-sphere and v and A are functions of u and r. The coordinate r is a 
generalization of the usual area radius. The coordinate u is constant on outgoing radial null geodesies. Here, we also 
require that on the world line of the central observer, u is equal to the proper time of that observer. This gives rise 
to the condition = at r = 0. Define Eab = Gab — STrTah where Gab is the Einstein tensor, so the Einstein-scalar 
equations are just Eab — 0. Note that as a consequence of the Bianchi identities, must satisfy the wave equation. 

In 4 spacetime dimensions, Christodoulou |^ showed how to write the spherically symmetric Einstein-scalar equa- 
tions as an integro-differential equation for the scalar field. We now generalize the method of reference j|] to the 
n-dimensional case. Define the null vectors l"^ and n° by 



or 



(4) 



d_ 

du 

Then some straightforward but tedious calculations give 

a,6 _ »^-2 



Gabl°'l 



2A 



dr 



dr 



(5) 



(6) 



Gabion'' 



One can show that the equations 



n - 2 



n-3 



,-2A 



(z.- A) 



Eabl^l'^Q , Eabl'^n'^^Q , V,V"(^ = 



(7) 



(8) 



plus the condition of regularity at the origin imply Eab — everywhere. Thus we need to impose only those three 
equations. 

Using equations (0,0) and (0), we find that the first two of equations (0) yield 



"-2 d ,^ , ^ {d(t> 
r or \ or 



(9) 



2 - n 
2r 



,-2A 



ri — 3 



1 - e 



-2X\ 



= 



(10) 



Define g = e'^^^ and g = e'^ ^. Then the solution of equations (O) and (Eoh is 
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g = exp 



(11) 



f"-%(f)df . (12) 







The wave equation for (j) in this metric is 

„ 73 - 2 d4> d 



--'^"tZ -""'ff^)=0 . (13) 



dudr r du dr \ dr 



In principle, equations ( |ll[]l3| ) could be used for a numerical treatment of the Einstein-scalar equations. Given 
on an initial light cone, equations (^Tj) and (^2|) could be used to find the metric on that light cone. Equation (|lJ 
could then be used to evolve to a nearby light cone. However, in 4 spacetime dimensions, Christodoulou ||^ finds 
that a change of variables gives rise to a nicer form of the wave equation. Define the operator D hy D = e'^n°'\7a- 
Then Z? is a derivative along ingoing light rays and we have 

Now in 4 spacetime dimensions define h = {d / dr){r(j)). Then the wave equation ( [T^ ) becomes 

Dh=^^{g - g){h -h) (15) 

where h = J^hdr. For the purposes of numerical simulations, the nice property of equation ( |l5| ) is that it involves 
a derivative only in the direction of ingoing light rays. Another nice feature of equation (|l^) is that in the absence 
of gravity {g — g — 1), the right hand side vanishes. This is related to the fact that in general, the right hand side 
vanishes like as r — > 0. (The metic is smooth, and therefore in a neighborhood of the origin it behaves like a flat 
metric to some order in r). This property is important for the following reason: spherically symmetric metrics have 
coordinate singularities at the origin. These coordinate singularities are reflected in the appearance of inverse powers 
of r in equation (^3|), which can lead to instabilities or inaccuracies in a numerical simulation. 

These considerations suggest that instead of a numerical simulation of equation (jlj) in n dimensions, we should 
instead search for a new variable h that satisfies Dh = in the absence of gravity. Note that a quantity that satisfies 
Dh = is constant along ingoing light rays and therefore satisfies Huygens' principle. However, it is well known 
that solutions of the wave equation in flat spacetime satisfy Huygens' principle only in even spacetime dimensions. |^ 
Therefore, we should only expect to find an appropriate new variable in even spacetime dimensions. Let n be even 
and define m = {n — 2)/2. For any (f) define 

- r^-fl-r (r^-V) . (16) 



(2m -1)! 

(Here, the numerical factor is chosen for later convenience). If (/) is a solution of the flat space wave equation (equation 
( |l3| ) with g = 1), then h satisfies Dh — 0. That is, in the absence of gravity, a solution of the wave equation (j) gives 
rise to an h that is constant along ingoing light rays. One can demonstrate this property of h by using equation ( [T^ ) 
(with g = 1) and mathematical induction on m. 

The advantage of using h as the basic variable is that it tends to make computer simulations more stable and 
accurate, especially near the origin. The disadvantages are (1) the method only works in even dimensions, (2) 
a seperate computer code must be written for each value of n and (3) as n gets larger, the equations become 
more complicated. We expect that some numerical technique can be used to evolve the Einstein-scalar equations 
in n dimensions with a single code, with n as a free parameter, stably and with enough accuracy to treat critical 
gravitational collapse. However, we have been unable to devise such a technique. Therefore, in this work we use h 
as our basic variable. Due to the increase in complication with increasing n, we treat only the case of 6 spacetime 
dimensions. 

We now specialize to the case of 6 dimensions. The variable h is given by 
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We must now use equation (11) to express g in terms of h and integrals involving h. Define the quantities s and s by 

2 



f h(f) df 



(18) 



s(f) df 



IStt / — (s - sf 





Then, using equations (|17|-|19D and (|ll| ) we find 

5 = exp 

Equation ( p^ specialized to 6 spacetimc dimensions becomes 

g = — g{f) df 

^ Jo 



(19) 



(20) 



(21) 



Using equations (17-2l|) in equation (|13| ) we find that the wave equation for (/) in 6 spacetimc dimensions becomes 

(22) 

Equations ([l7j-|22|) are the full set of equations that are evolved in 6 spacetimc dimensions. 



Dh =- {g - g){h + S - 2s) + ^ {S - sf 

r 2 r 



IV. NUMERICAL METHOD 

The numerical method used is that of Garfinkle [^. This method is based on numerical work of Goldwirth and 
Piran , which is in turn based on the analytical work of Christodoulou . The spatial grid is a set of points on an 
outgoing null cone, and each spatial grid point is evolved along an ingoing light ray. Given /i on a light cone, the code 
performs several integrals in turn: equations ( p^|^l| ). Near the origin, these integrals are approximated by a Taylor 
series and evaluated using the slope of ft, at r = 0. At all other points, the integrals are evaluated using Simpson's 
rule. 

For each time step the quantity h is evolved using equation (|2^ ) and the quantity r is evolved using 



Dr = - 



(23) 



These are essentially a set of uncoupled ODEs, one for each grid point. In the scheme of |^ as the evolution precedes 
grid points that pass through r = are lost. When half of the grid points are lost, they are put back, interpolated 
between the remaining grid points. The critical solution is found by a search of a 1-parameter space of evolved data to 
find the boundary between those data that form black holes and those that do not. The outermost gridpoint is chosen 
to be the light ray that hits the singularity of the critical solution. This choice maintains resolution throughout the 
evolution. 



V. RESULTS 

All runs were done with 300 spatial grid points. The code was run in quadruple precision on Dec alpha workstations 
and in double precision on a Cray YMP8. The initial data for the scalar field was chosen to be of the form 



0(0, r) = p exp — (r — ro) /cr 



(24) 



Here, p is our parameter, and rg and a are constants. The critical solution was found to have discrete self-similarity. 
Let u* be the value of u at the origin at the singularity of the critical solution. Define the coordinates T = — ln(M* — u) 



and R = 



Then discrete self similarity means that h, considered as a function of T and R, is periodic in T. Figure 



1 shows h plotted as a function of T and R. Note that after some initial transient behavior, h becomes periodic. To 
examine the periodicity more closely, in Figure 2 we plot h vs R bX two different times where the minimum of h occurs 
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at r = 0. Note that the curves agree, demonstrating periodicity of h. The period of /i is A = 3.03. Thus, in each 
period a typical length of the system shrinks by a factor of « 21. For comparison, recall that for n = 4, A = 3.445. 

To compute the scaling exponent 7 we use the method of reference [^. We evolve data for a range of parameters 
below the threshold of black hole formation. For each evolution we find the maximum of the absolute value of the 
scalar curvature on the world line of the central observer Rmax- Plotting lni?„iax vs ln(p * —p) the result is a straight 
line with a periodic wiggle, where the slope of the line is — 27. We use 50 values of p equally spaced in ln(p * —p)- 
Figure 3 is a plot of lni?max vs ln{p * —p) along with the best straight line fit to the points. The slope allows us to 
find 7. The result is 7 = 0.424, so Mbh oc {p — p*)^'^^. Recall that for n = 4, 7 = 0.374. Note that 7 can also be 
found by evolving data above the threshold of black hole formation and plotting In Mbh vs \n{p~p*). We have done 
this and the result is a straight line with a periodic wiggle where the slope of the line is w 1.27. However, our method 
gives a more accurate treatment of subcritical collapse than of supercritical collapse. For that reason, we have used 
subcritical collapse to calculate 7. 

Comparison of our results for 6 dimensions with Choptuik's results for 4 dimensions 0| seems to indicate that A is 
a decreasing function of the dimensionality n of spacetime. The quantity 7 seems to be an increasing function of n. 
It would be interesting to obtain more information about these functions by studying critical collapse for other values 
of n. 
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FIG. 1. h plotted vs. the rescaled coordinates R and T. 
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FIG. 2. h plotted vs. R at two times when the minimum of ft is at r = 0. Note that the two curves agree. 
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FIG. 3. In -Rmax plotted vs. \n{p * —p) along with the best straight line fit. The curve is a straight line with a periodic wiggle 
and the slope of the line is —27. 
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